#############################################################################################
## Replication code for:                                                                   ##
## A Male Hostility Spiral? Polarized Communication among Political Elites on Social Media ##
## Step 7.                                                                                 ##
#############################################################################################

library(reshape2)
library(ggeffects)
library(igraph)
library(ggraph)
library(patchwork)
library(progress)
library(tidyverse)
options(scipen=999)

data <- read.csv("data1.csv")

data$predicted_targeted_sentiment <- recode(data$predicted_targeted_sentiment,
                                            "Negative" = -1,
                                            "Neutral" = 0,
                                            "Positive" = 1)

data$inparty <- as.numeric(data$party_id_receiver==data$party_id_sender)
data$sender_gender <- as.factor(data$sender_gender)
data$receiver_gender <- as.factor(data$receiver_gender)

ols_model <- lm(predicted_targeted_sentiment~inparty*sender_gender*receiver_gender,
                data=data)

one_iteration <- function(n_pol=30,
                           num_step=20,
                           prop_men=0.5,
                           party_prop=c(0.5,0.5),
                           sentiment_lag1=0.2,
                           decay=1,
                           model=ols_model){
  
  num_parties = length(party_prop)
  genders <- c(rep("male", round(prop_men*n_pol)), rep("female", round((1-prop_men)*n_pol)))
  genders = sample(genders)
  
  parties = c(rep("partyA",party_prop[1]*n_pol),rep("partyB",party_prop[2]*n_pol))
  parties <- sample(parties)
  
  A = data.frame(id=rep(NA,n_pol),
                 party=rep(NA,n_pol),
                 gender=rep(NA,n_pol))
  
  A$id = 1:n_pol
  A$party <- parties
  A$gender <- genders
  
  S <- matrix(0,nrow=n_pol,ncol=n_pol)
  rownames(S) <- 1:n_pol
  colnames(S) <- 1:n_pol
  
  
  for (i in 1:n_pol){
    for (j in 1:n_pol){
      observation <- data.frame(sender_gender = A$gender[i],
                                receiver_gender = A$gender[j],
                                inparty = as.numeric(A$party[i]==A$party[j]))
      S[i,j] <- predict(model,observation)+rnorm(1,0,sd(model$residuals))

    }
  }
  
  diag(S) <- NA
  S <- tanh(S)
  S_iter <- list()
  
  for (t in 1:num_step){
    S_iter[[t]] <- S
    
    for (i in 1:n_pol){
      for (j in 1:n_pol){
        S[i,j] <- atanh(S[i,j])
        
        prev_int <- atanh(S_iter[[t]][j,i])

        lr <- sentiment_lag1*decay^t
        S[i,j] <- tanh(S[i,j]+lr*prev_int)
      }
    }
    
  }
  return(list(S_iter,A))
  
}

sim_1 <- function(iteration){
  S_iter <- iteration[[1]]
  A <- data.frame(iteration[[2]])
  dyadic_df_list <- list()
  
  for (i in seq_along(S_iter)) {
    sentiments <- S_iter[[i]]
    
    sentiments_long <- as.data.frame(sentiments) %>%
      rownames_to_column(var = "Sender_ID") %>%
      pivot_longer(cols = -Sender_ID, names_to = "Receiver_ID", values_to = "Sentiment")
    sentiments_long$Sender_ID <- as.numeric(sentiments_long$Sender_ID)
    sentiments_long$Receiver_ID <- as.numeric(sentiments_long$Receiver_ID)
    
    dyadic_df <- sentiments_long %>%
      left_join(A, by = c("Sender_ID" = "id")) %>%
      rename(Sender_Party = party, Sender_Gender = gender) %>%
      left_join(A, by = c("Receiver_ID" = "id")) %>%
      rename(Receiver_Party = party, Receiver_Gender = gender) %>%
      select(Sender_Party, Sender_Gender, Receiver_Party, Receiver_Gender, Sentiment)
    
    dyadic_df_list[[i]] <- dyadic_df
  }
  
  
  average_sentiments_df <- data.frame(
    Time_Step = integer(),
    Inparty = numeric(),
    Outparty = numeric()
  )
  
  for (i in seq_along(dyadic_df_list)) {
    dyadic_df <- dyadic_df_list[[i]]
    
    avg_same_party <- dyadic_df %>%
      filter(Sender_Party == Receiver_Party) %>%
      summarise(Avg_Sentiment = mean(Sentiment, na.rm = TRUE)) %>%
      pull(Avg_Sentiment)
    
    avg_different_party <- dyadic_df %>%
      filter(Sender_Party != Receiver_Party) %>%
      summarise(Avg_Sentiment = mean(Sentiment, na.rm = TRUE)) %>%
      pull(Avg_Sentiment)
    
    average_sentiments_df <- rbind(average_sentiments_df, data.frame(
      Time_Step = i,
      Inparty = avg_same_party,
      Outparty = avg_different_party
    ))
  }
  
  result <- melt(average_sentiments_df,id.vars="Time_Step")
  
  return(result)
}

sim_2 <- function(iteration){
  S_iter <- iteration[[1]]
  A <- data.frame(iteration[[2]])
  dyadic_df_list <- list()
  
  for (i in seq_along(S_iter)) {
    sentiments <- S_iter[[i]]
    
    sentiments_long <- as.data.frame(sentiments) %>%
      rownames_to_column(var = "Sender_ID") %>%
      pivot_longer(cols = -Sender_ID, names_to = "Receiver_ID", values_to = "Sentiment")
    sentiments_long$Sender_ID <- as.numeric(sentiments_long$Sender_ID)
    sentiments_long$Receiver_ID <- as.numeric(sentiments_long$Receiver_ID)
    
    dyadic_df <- sentiments_long %>%
      left_join(A, by = c("Sender_ID" = "id")) %>%
      rename(sender_party = party, sender_gender = gender) %>%
      left_join(A, by = c("Receiver_ID" = "id")) %>%
      rename(receiver_party = party, receiver_gender = gender) %>%
      select(sender_party, sender_gender, receiver_party, receiver_gender, Sentiment)
    
    dyadic_df_list[[i]] <- dyadic_df
  }
  
  new_data <- expand.grid(
    inparty = c(0, 1),
    receiver_gender = c("male", "female"),
    sender_gender = c("male", "female")
  )

  predicted_effects_joint <- NA
  for (i in 1:length(dyadic_df_list)) {
    dyadic_df <- dyadic_df_list[[i]]
    dyadic_df$inparty <- as.numeric(dyadic_df$sender_party == dyadic_df$receiver_party)
    dyadic_df$sender_gender <- factor(dyadic_df$sender_gender,levels=c("male","female"))
    dyadic_df$receiver_gender <- factor(dyadic_df$receiver_gender,levels=c("male","female"))
    
    model <- lm(Sentiment ~ inparty*sender_gender*receiver_gender, data = dyadic_df)
    
    predicted_effects <- data.frame(
      new_data,
      predicted_value = predict(model, newdata = new_data)
    )
    

    colnames(predicted_effects) <- c("Inparty","Receiver Gender","Sender Gender","estimate")
    predicted_effects$Inparty <- factor(predicted_effects$Inparty)
    levels(predicted_effects$Inparty) <- c("Outparty","Inparty")
    predicted_effects$`Sender Gender` <- factor(predicted_effects$`Sender Gender`,levels=c("female","male"))
    levels(predicted_effects$`Sender Gender`) <- c("Sender Female","Sender Male")
    predicted_effects$`Receiver Gender` <- factor(predicted_effects$`Receiver Gender`)
    levels(predicted_effects$`Receiver Gender`) <- c("Receiver Male","Receiver Female")
    predicted_effects$iter <- i
    predicted_effects_joint <- rbind(predicted_effects_joint,predicted_effects)
  }
  predicted_effects_joint <- predicted_effects_joint[-1,]
  return(predicted_effects_joint)
}


sim_3 <- function(iteration){
  S <- iteration[[1]][[length(iteration)]]
  A <- iteration[[2]]
  A$party <- factor(A$party)
  levels(A$party) <- c("A","B")
  A$gender <- factor(A$gender,levels=c("female","male"))
  
  sentiments_long <- as.data.frame(S) %>%
    rownames_to_column(var = "Sender_ID") %>%
    pivot_longer(cols = -Sender_ID, names_to = "Receiver_ID", values_to = "Sentiment")
  sentiments_long$Sender_ID <- as.numeric(sentiments_long$Sender_ID)
  sentiments_long$Receiver_ID <- as.numeric(sentiments_long$Receiver_ID)
  
  dyadic_df <- sentiments_long %>%
    left_join(A, by = c("Sender_ID" = "id")) %>%
    rename(Sender_Party = party, Sender_Gender = gender) %>%
    left_join(A, by = c("Receiver_ID" = "id")) %>%
    rename(Receiver_Party = party, Receiver_Gender = gender)
  dyadic_df <- dyadic_df[!is.na(dyadic_df$Sentiment),]
  edge_list <- data.frame(from = dyadic_df$Sender_ID, to = dyadic_df$Receiver_ID, weight = dyadic_df$Sentiment+2)
  g <- graph_from_data_frame(edge_list, directed = TRUE)
  
  
  for (i in 1:nrow(A)){
    V(g)[i]$Party <- A[i,]$party
    V(g)[i]$Gender <- A[i,]$gender
  }
  
    
  return(g)
}

sim_1_ci <- function(num_iter,male_prop){
  pb <- progress_bar$new(
    format = "[:bar] :percent ETA: :eta",
    total = num_iter,
    clear = FALSE,
    width = 60
  )
  temp <- vector("list", num_iter)
  
  for (i in 1:num_iter){
    pb$tick()
    temp[[i]] <- sim_1(one_iteration(prop_men=male_prop))
  }
  
  temp_df <- do.call(rbind, temp)
  
  output <- temp_df %>%
    group_by(Time_Step, variable) %>%
    summarise(
      mean_value = mean(value),
      ci.high=mean(value)+1.96*sd(value),
      ci.low=mean(value)-1.96*sd(value),
      .groups = "drop"
    )
  output <- as.data.frame(output)
  return(output)
}

graph_data_combined <- NA
male_props <- c(0,0.25,0.5,0.75,1)
for (i in 1:5){
  graph_data <- sim_1_ci(num_iter = 300,
                         male_prop=male_props[i])
  graph_data$prop_men <- male_props[i]
  graph_data_combined <- rbind(graph_data_combined,graph_data)
  rm(graph_data)
}
graph_data_combined <- graph_data_combined[-1,]

png(file="tables_and_figures/fg5.png",
   width=16, height=10, units="cm", res=600)
ggplot(graph_data_combined,aes(Time_Step,mean_value,col=variable,group=variable))+
  geom_line()+
  geom_line(aes(x=Time_Step,y=ci.high,col=variable),lty="dashed")+
  geom_line(aes(x=Time_Step,y=ci.low,col=variable),lty="dashed")+
  facet_wrap(~prop_men,nrow=1)+
  theme_bw()+
  xlab("Iterations")+
  ylab("Average Sentiment")+
  theme(legend.position="bottom",
        legend.title=element_blank())
dev.off()

male_props <- c(0,1)
network_graphs <- list()
for (i in 1:length(male_props)){
  res <- one_iteration(prop_men=male_props[i],decay=1,n_pol=40,num_step=30)
  g <- sim_3(res)
  
  layout <- layout_with_fr(g, 
                           weights = abs(E(g)$weight) * 20,
                           niter = 3000,
                           grid = "nogrid")
  
  network_graphs[[i]] <- ggraph(g, layout = layout) +
    geom_edge_link(aes(color = weight, alpha = abs(weight)/7-0.3), width = 0.4, show.legend = FALSE) +
    geom_node_point(aes(color = Party,shape=Gender), size = 3,show.legend=F) +
    scale_edge_color_gradient(low = "red", high = "green") +
    scale_shape_manual(values = c(17, 16), labels = c("female", "male"), drop = FALSE) +
    theme_bw() +
    scale_color_manual(values = c("red", "blue")) + 
    theme(legend.position = "bottom") +
    ylab("") +
    xlab("")+ 
    theme(axis.text.x=element_blank(), 
          axis.ticks.x=element_blank(), 
          axis.text.y=element_blank(), 
          axis.ticks.y=element_blank())+
    theme(plot.title = element_text(hjust = 0.5))
}

graph_plot <- wrap_plots(network_graphs, ncol =length(male_props)) +
  plot_layout(guides = 'collect')& theme(legend.position = "bottom")

png(file="tables_and_figures/fg6.png",
    width=20, height=8, units="cm", res=600)
graph_plot
dev.off()


sim_2_ci <- function(num_iter){
  pb <- progress_bar$new(
    format = "[:bar] :percent ETA: :eta",
    total = num_iter,
    clear = FALSE,
    width = 60
  )
  temp <- vector("list", num_iter)
  
  for (i in 1:num_iter){
    pb$tick()
    temp[[i]] <- sim_2(one_iteration(prop_men=0.5,n_pol=30,decay=1,num_step=20))
  }
  
  temp_df <- do.call(rbind, temp)
  

  output <- temp_df %>%
    group_by(Inparty, `Sender Gender`, `Receiver Gender`,iter) %>%
    summarise(
      mean_estimate = mean(estimate),
      ci.high=mean(estimate)+1.96*sd(estimate),
      ci.low=mean(estimate)-1.96*sd(estimate),
      .groups = "drop"
    )
  
  return(output)
}
sim_2_graph <- sim_2_ci(num_iter=300)

png(file="tables_and_figures/fg7.png",
    width=12, height=12, units="cm", res=600)
ggplot(sim_2_graph, aes(x=iter,y = mean_estimate,col=Inparty,group=Inparty)) +
  geom_line()+
  geom_line(aes(x=iter,y=ci.high,col=Inparty),lty="dashed")+
  geom_line(aes(x=iter,y=ci.low,col=Inparty),lty="dashed")+
  facet_grid(`Sender Gender`~`Receiver Gender`) +
  theme_bw()+
  xlab("")+
  ylab("Predicted Sentiment")+
  theme(legend.position="bottom")
dev.off()

